Double occupancies in confined attractive fermions on optical lattices 
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We perform a numerical study of a one-dimensional Fermion-Hubbard model in harmonic 
traps within the Thomas-Fermi approximation based on the exact Bethe-ansatz solution. The 
p — U/t phase diagram is shown for the systems of attractive interactions (p is the charac- 
teristic density and U/t the interaction strength scaled in units of the hopping parameter.). 
We study the double occupancy, the local central density and their derivatives. Their roles 
are discussed in details in detecting the composite phases induced by the trapping potential. 
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The recent remarkable experimental progress in cool- 
ing and manipulating ultracold atomic gases in optical 
lattices provides us an alternative way to investigate the 
many-body effects in condensed matter system [1|. For 
the optical lattices loaded with atomic gases, the on-site 
interaction between different species can be tuned to be 
very strong by means of a technique named Feshbach 
resonance Q or by increasing the lattice depth to de- 
crease the hopping between neighbor sites. Thus the 
outstanding controllability offers the possibility to use 
these systems to "quantum simulate" quantum many- 
body physics in a well-controlled way. An example is the 
experimentally realized Tonks- Girardeau gas in a one- 
dimensional (ID) optical lattice Q. 

In the experiments, the trapping potential used to 
confine atoms induces inhomogeneity and complexity in 
characterizing the quantum phases. Though a few tech- 
niques exist in experiments observing the new quantum 
phases and probing various properties in ultracold gases 
in optical lattices, such as the time-of-fiight imaging of 
the momentum distribution, noise correlations and Bragg 
spectroscopy there are difficulties in observing the 
composite phases induced by the inherent trapping po- 
tential, for example, the incompressible bulk Mott- and 
band-insulating phases surrounded by the compressible 
fluid phases at the edges in the ID Fermi-Hubbard model 
under harmonic confinement. Scarola et al. proposed to 
use the core compressibility by filtering out the edge ef- 
fects, which offers a direct probe of incompressible phases 
independent of inhomogeneity Q. Duan et al. suggested 
to use the Fourier sampling of time-of-flight images to re- 
veal new correlation functions [|| . Kollath et al. used the 
double occupancy (DO) induced by periodic lattice mod- 
ulations via the creation of molecules to identify the pair- 
binding energy and the spin-ordering in a Fermi gas in 
an optical lattice @ , which was recently measured in the 
experiments by molecule formation and used to identify 
an incompressible Mott-insulator phase of two hypcrfinc 
states of strongly repulsive fermionic 40 K atoms Q- The 



compressibility of the quantum gas in the trap was used 
recently in identifying the incompressible Mott-insulating 
phase at strong interactions by independent control of 
external confinement and lattice depth [8|. Here in this 
paper we discuss how to characterize the different com- 
posite phases of ID confined fermions on optical lattices, 
by investigating observables like the double occupation 
sites, the local central density and their derivatives. 

We consider a two-component Fermi gas in a tube with 
Nf atoms and N s lattice sites, which can be realized 
by a strong confinement in the transverse directions Q 
with an additional periodic potential applied along the 
tube and can be described by a one-band inhomogeneous 
Fermi-Hubbard model [Toj . 
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where the spin degrees of freedom a =\, 4- ar e pseudospin- 
1/2 labels for two internal hyperfine states. Ci a (y icr ) are 
fermionic operators annihilating (creating) particles with 
spin cr in a Wannier state at site i. fii = hi a = 
J2 a Hosier is the total site occupation operator, t is the 
tunneling between nearest neighbors, U is the strength 
of the on-site attractive interaction, Di = hi-fhii is the 
double occupancy operator at site i, and Vi describes the 
strength of the trapping potential. 

The homogeneous ID Fermi-Hubbard model belongs 
to the universality class of Luttinger liquids. At zero 
temperature, the properties of this model in the ther- 
modynamic limit (N a ,N s — !• oo, but with finite N a /N s ) 
are determined by the fillings n a — N a /N s and by the 
dimensionless coupling constant u = U/t. According to 
Lieb and Wu the ground state properties for differ- 
ent fillings in the thermodynamic limit are described by 
the coupled integral equations (details see Refs. [H, [H). 

The inhomogeneity of ID Fermi-Hubbard model 
caused by the boundary conditions or the external po- 
tential, normally invalidates a reliable analytical method 
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FIG. 1: (Color online) Phase diagram for the one-dimensional 
inhomogeneous fermion-Hubbard model as a function of at- 
tractive interaction strength u = U/t and characteristic par- 
ticle density p = Nf^/Vz/i. For fixed u, the shape of the 
density profile is completely determined by p. The band in- 
sulator (BI) phase is determined with the local occupation 
in the trap center satisfying \m — 2| < 0.005. Below the BI 
phase, it is the Luther-Emery (LE) liquid phase. The bound- 
ary between BI and LE phases is plotted as open squares (□) 
with solid lines as a guide for the eye. The solid circles (•) 
represent the phase boundary determined from the derivative 
of the double occupancy, and the solid upper triangles (a) 
represent the one determined from the derivative of the local 
central density (See details in the text and in Fig. All of 
them can be used as identifying the onset of the BI phase. 




that usually used in the homogeneous system. Many nu- 
merically accurate schemes, such as the density-matrix 
renormalization group EH, E3, EH , quantum- Monto-carlo 
(QMC) EMU, exact diagnolization El, US H3 and 
density-functional theory based on the exact Bethe- 
ansatz solution E2], are used in obtaining the phase di- 
agram and the collective oscillations of the atomic mass 
density. The Thomas-Fermi approximation based on the 
exact Bethe-ansatz solution has also been successfully 
used in characterizing the phase diagrams [22j . The rich 
phases induced by the trapping potential or the phase 
transitions driven by the system parameters can be iden- 
tified by the site occupation, the variance of the local 
density, the local compressibility [l9j |. or the thermody- 
namic stiffness (i.e., the inverse of the global compress- 
ibility) El- 

Here we study the inhomogeneous Fermi-Hubbard 
model described by Eq. ([T]), which allows for the spa- 
tial coexistence of different quantum states [18j, |22| . For 
attractive interactions, U < 0, there are two different 
composite phases, well classified in the phase diagram 
defined b y a s caled dimensionless characteristic density 
P = N fv / V 2 /t and u = U/t 18]: the Luther-Emery (LE) 
liquid phase characterized by fermion pairs and atomic 



FIG. 2: (Color online) The ground-state concentration of dou- 
bly occupied sites in the homogeneous system £) hom ( a ) as a 
function of interaction strength u for various n and (b) as a 
function of particle density n for various interacting strengths 
u. For < n < 2, D hom is in the range < D hom < I. For 
fixed n, £) hom increases monotonically with increasing |m| and 
for fixed u it increases monotonically with increasing n. At 
U = 0, D houl — n 2 /4 and at strongly attractive coupling 
limit D houl — n/2. The attractive interaction favors the for- 
mation of pairs between opposite spins. In the inset of (b), 
the derivative of the double occupancy with respect to the 
particle density n is plotted, which is used to calculate the 
double occupancy for the inhomogeneous system in the spirit 
of local density approximation. 



density wave, and band insulator (BI) phase with a 
plateau of n.; = 2 in the trap center of tightly bound spin- 
singlet dimers surrounded by Luther-Emery layers I12|. 

All these phases have been analyzed in the p — u |19l . 
I23L u — u (fi is the chemical potential) J24j, and n — 
v [22| (no is the central density and v = 2^J Nf /it is the 
characteristic filling factor) phase diagrams. 

The harmonic trapping potential in Eq. ([1]) induces 
a position dependence of the chemical potentials which 
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FIG. 3: (Color online) (a) The ground-state concentration of 
doubly occupied sites D as a function of characteristic density 
p in the system under the presence of trapping potential for 
various interacting strength it. (b) The local central density 
no as a function of characteristic density p for various inter- 
acting strength u. For comparison, the double occupancy and 
the local central density for the noninteracting case (u = 0) 
are also shown in (a) and (b), respectively. 



can be evaluated within the local density approximation 
based on the Bethe-ansatz equations. This amounts to 
determining the global chemical potential p of the system 
from the local equilibrium condition, 

M = /ihom(n,u)| n ^ n . +V 2 {i- N s /2) 2 , (2) 

while imposing the normalization condition rii = Nf. 
Here /ihom(^)W) is the chemical potential of the homoge- 
neous system calculated from the Bethe-ansatz coupled 
equations [r2T l25j. 

In Fig. 1, we show the above-mentioned composite 
phases through the p — u phase diagram, which is insen- 
sitive to the individual values of Nf and V2, that is, for 
fixed u, the shape of the density profile in this system is 
completely determined by p [l8(. Below the BI phase, it 
is the LE liquid phase. We would like to mention that the 
similar phase diagram has been obtained by Rigol et al. 
for the repulsive interaction [l8[ and Heiselberg but with 
a fitting function for Hbom(n,u) [23J. Here we solve the 
local equilibrium condition with the exact numerical data 
for /Xhom(^)W) calculated from the coupled Bethe-ansatz 
equations. We conclude that the phase diagram based 
on the Thomas-Fermi approximation gives qualitatively 
correct results when comparing to QMC. 

To detect the composite phases induced by the trap- 
ping potential, we investigate firstly the concentration of 
the doubly occupied sites, which is defined as the deriva- 
tive of the ground state energy with respect to U 26] , 

Nf ^ 1 N f dU 

i— 1 J 

Here E is the ground state energy per site of the sys- 
tem. The DO in the homogeneous system, D hom , can be 
calculated exactly from the Bethe-ansatz coupled equa- 
tions [l3j], which is shown in Fig. 2. Notice, at U = 0, we 
have D hom = n 2 /A and at strongly attractive coupling 
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FIG. 4: (Color online) (a) The derivative of the ground-state 
concentration of doubly occupied sites D as a function of char- 
acteristic density p in the system under the presence of trap- 
ping potential. In the inset, the local inverse compressibility 
n^ 1 (in units of t) is shown as a function of i with the sys- 
tem parameters determined from the critical characteristic 
density on the peaks, (b) The derivative of the local central 
density no as a function of characteristic density p for differ- 
ent attractive interactions. For comparison, the derivatives 
of the double occupancy and the local central density for the 
noninteracting case (u = 0) are also shown in (a) and (b), 
respectively. 



limit £) hom = 71/2. In the presence of trapping poten- 
tial, we calculate the DO in the spirit of local density 
approximation, D = -j^ ^\ D hora (n)\ n ^ ni . From Fig. 2, 
it becomes obvious that the attractive interaction favors 
the formation of pairs between the two spin species while 
the repulsive interaction suppresses them. 

Now we show in Fig. 3 (a) the DO and in Fig. 3 
(b) the local central density no for the inhomogeneous 
system described in Eq. (1) as a function of the charac- 
teristic density p. For large p, the DO and the central 
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density saturate at D w 1 and no w 2, respectively, cor- 
responding to the BI regime, where an incompressible 
phase of band-insulating type is realized in the core of 
the trap. The more attractive of the interaction is, the 
sooner for D to reach the value of 1 and for no to reach 
2. These manifest that the attractive interaction favors 
pairing, and thus increases DO, different from the case 
of repulsive interaction where it suppresses DO and sup- 
ports Mott-insulating phase. 

In the following, we show the derivative of the DO in 
Fig. 4 (a) and the derivative of the local central density 
in Fig. 4 (b), respectively, as probes for the compos- 
ite quantum phase transitions and especially as signs for 
band-insulating phases. The derivatives of the DO and 
the local central density are performed numerically us- 
ing the central finite difference method. We find that 
the singularities in these two quantities provide a further 
support and clear signatures of the composite quantum 
phase transitions accompanying Fig. 3. In Fig. 4 (a), 
the phase transition from LE to BI is manifested by an 
obvious peak. The onset of the peak gives a critical value 
of the characteristic density, which determines the crit- 
ical particle number for a given Vi/t. Based on these 
parameters, we plot in the inset of Fig. 4 (a) the local 
inverse compressibility defined by n^ 1 — d 2 Eo/dn 2 \ n ^ ni 
as a function of site position. The incompressible nature 
of the BI phase is manifested by the disappearance of 



k~ , which is compatible with the onset of the peak in 
dD/dp. In Fig. 4 (b), the vanishing of the derivative 
of the local central density gives very direct signatures 
of the incompressible nature of BI in the center of the 
trap. The phase boundaries of the onset of the BI phase 
are plotted in Fig. [TJ which are determined by the ap- 
pearance of the peak in the derivative of the DO and 
the vanishing of the derivative of the local central den- 
sity. We find that the phase boundary determined from 
the DO slightly overestimates the onset of the BI phase, 
while the central density slightly underestimates it. 

In conclusion, motivated by earlier theoretical work 
and by ongoing experimental efforts, we present a fully 
numerical study for the ground state properties of a 
one-dimensional attractive Fermion-Hubbard model in 
harmonic traps within the Thomas-Fermi approximation 
based on the exact Bethe-ansatz solution. We investi- 
gate the double occupancy, the local central density and 
their derivatives. We found they are good indicators of 
the formation of an incompressible band-insulating phase 
in the center of the trap, which can be detected in the 
near-future in the experiments in the ID cold fermionic 
atomic gases. 
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